% ESTIMA_LOOP:
% Function can also be called with one giant function handle, may be
% slower.
% estfunc = @(x)lmj_postmin(x, parstrval, parposest, ...
%                  partrmat_est, funcmod, Y, trainvec, prpar, prss, ...
%                  solveopt, addsol, ssposest, flag_transform);
%
% NUMPAR            Number of total parameters (including calibrated)

% 1. Possibly start a DIARYFILE
% =========================================================================
if flag_diary==1
    diaryfname=['log ',strdate(1)];
    % Create a Diary File of the Optimization
    delete(diaryfname);
    diary(diaryfname);
    disp('__________________________________________________');
    datestr(now);
end

% 2.Organize Outputs
% =========================================================================
logprior_mat(:,2)= -1e15*ones(nstrvals,1);
loglikel_mat(:,2)= -1e15*ones(nstrvals,1);
logpostd_mat(:,2)= -1e15*ones(nstrvals,1);
% PARMODEMAT     Matrix of output vectors
% OPTIMAT        Matrix with info on optimization
%                First columns number of iterations
%                Second column return flags
%                Third time in minutes
parmodemat = zeros(numpar, nstrvals);
optimat    = zeros(3, nstrvals);

% 2.
% Begin estimation Loop
% Will work with TRY and CATCH later on but not yet for debugging

% 3. Start a cell that can save output as we go along
% TAB_LOOP
% =========================================================================
tab_loop=emptycell(numpar+6,nstrvals+1);
tab_loop(1:numpar,1)=parnames;
tab_loop(numpar+1:end,1)={'logPost';'logLikel';'Number';'Iteration count';'exitflag';'minutes'};
cd(outpath);save tab_loop;cd(cucd);


% 4. Assign default values to ESTIMOPT if they do not exist
% =========================================================================
if ~exist('estimopt', 'var'); estimopt = []; end;
[junk,estimopt]=ch_field(estimopt,'TolFun',1e-5);
[junk,estimopt]=ch_field(estimopt,'MaxFuneval',200000);
[junk,estimopt]=ch_field(estimopt,'MaxIter',400);
[junk,estimopt]=ch_field(estimopt,'Display','on');
if case_estima == 4; % when using Simann
    % pass in special inputs
    [junk,estimopt]=ch_field(estimopt,'sa_t',1e7);
    [junk,estimopt]=ch_field(estimopt,'sa_rt',0.9);
    [junk,estimopt]=ch_field(estimopt,'sa_nt',5);
    [junk,estimopt]=ch_field(estimopt,'sa_ns',20);
    [junk,estimopt]=ch_field(estimopt,'rseed',999);
    [junk,estimopt]=ch_field(estimopt,'sa_neps', 4);
    [junk,estimopt]=ch_field(estimopt,'sa_eps', 1e-5); % convergenece criteria
    [junk,estimopt]=ch_field(estimopt,'sa_maxeval', 12*1e6);
    [junk,estimopt]=ch_field(estimopt,'sa_maxiter',  3*1e6); %3*1e6);
end

% 5. Assign default values to a PATH in case connection breaks down
% ------------------------------------------------------------------------
if isunix==0
    backup_path=cr_dir('C:','matlabrepo',root,spec,subf);
else
    backup_path=outpath;
end
% =========================================================================
% 6.
% Preliminary assigments and truncations before loop begins
% Formerly known as INDMAT, this is the matrix of transformations
% PARVECMODE is an empty vector with the calibrated values in
% =========================================================================
partrmat_est          =partrmat(parposest, :);
parvecmode            =zeros(numpar, 1);
parvecmode(parposcal) =parcalval;
errcount              =0;

% Close all paralell jobs 
matlabpool close force
matlabpool(8);

parfor ii=1:nstrvals;
    
    disp('===================================================');
    disp('===================================================');
    disp(['Starting Candidate Value number : ',num2str(ii) ] );
    xzero=parstrvals(parposest,ii);
    flag_optimok=1;
    
    tic;

% ========================================================================
% Output that must be produced by each code
% -----------------------------------------
% XESTMODE:  Estimated value
% ESTIMITER: Number of Iterations
% EXITFLAG : How routine ended
% LPOSTMIN : Exit value of posterior (min)
% April 19 2010
% ========================================================================
try
    switch case_estima
        case 1
            
            % =================================================================
            % 1. CSMINWEL
            % IO: transform from MOD 2 MIN and viceversa
            % =================================================================
            xzeromin          = mod2min(xzero,partrmat_est);
            
            [lpostmin,xestmin,junk,H,estimiter,junk,exitflag]=...
                csminwel('lmj_postoptim_central',xzeromin, 0.25*eye(length(xzeromin)),[],estimopt.TolFun,estimopt.MaxIter, ...
                parvecmode,parposest,partrmat_est,funcmod,Y,trainvec,prpar,prss,solveopt,addsol,ssposest,1);
            
            xestmode  = min2mod(xestmin,partrmat_est);
            
        case 2 % FMINCON
            
            % =================================================================
            % 2. FMINCON
            % ESTIMA_FMINCON.m
            % IO No transformation needed
            % Bounds transformed inside code
            % =================================================================
            [lpostmin,xestmode,exitflag,fminconout,Hx]=estima_fmincon_central(xzero,parvecmode,parposest,funcmod,Y,trainvec,prpar,prss,solveopt,addsol,ssposest,estimopt);
            estimiter = fminconout.iterations;
            % Old way of calling FMINCON; use subfunction directly
            %             flag_transform = 0;
            %             [xestmode, junk, exitflag, fminconout] = fmincon(@(x)lmj_postoptim(x,parstrval,parposest, ...
            %                 partrmat_est, funcmod, Y, trainvec, prpar, prss, solveopt, addsol, ssposest, flag_transform), ...
            %                 x0mode,[],[],[],[],prpar_lb(parposest),prpar_ub(parposest),[], estimopt);
            %             estimiter = fminconout.iterations;
            % [xestmode, junk, exitflag, fminconout] = fmincon(estfunc, x0mode,[],[],[],[],prpar_lb(parposest),prpar_ub(parposest),[], estimopt);
        case 3      % Giorgio's SA
            % =================================================================
            % 3. Simple version of SA
            % IO: transform from MOD 2 MIN and viceversa
            % ================================================================
            xzeromin          = mod2min(xzero,partrmat_est);
            
            nest=4*round(estimopt.MaxIter/4);
            dispaj('Running SA for ',nest,' iterations');
            [fh,xestmin,estimiter,exitflag,fvec]=simulated_simp(@lmj_postoptim,xzeromin(:) ...
                ,2,1e5,nest,parvecmode,parposest,partrmat_est,funcmod,Y,trainvec,prpar,prss,solveopt,addsol,ssposest,1);
            
            xestmode  = min2mod(xestmin,partrmat_est);
            
        case 4
            % =================================================================
            % 4. SIMMAN
            % IO: No transformation needed, except bounds passed from
            %     outside
            % ================================================================
            
            [xestmode, exitflag, estimiter, fvec]=simann(@(x)lmj_postoptim(x, parvecmode, parposest, ...
                partrmat_est, funcmod, Y, trainvec, prpar, prss, solveopt, addsol, ssposest, 0), ...
                xzero, prpar_lb(parposest),prpar_ub(parposest), estimopt);
            
        case 5
            
            % =================================================================
            % 4. FMINSEARCH
            % ESTIMA_FMINCON.m
            % IO: No transformation needed, bounds truncated inside code
            %     outside
            % =============================================================
            [lpostmin,xestmode,exitflag,fminsearchout]=estima_fminsearch(xzero,parvecmode,parposest,partrmat_est,funcmod,Y,trainvec,prpar,prss,solveopt,addsol,ssposest,estimopt);
            flag_transform = 1;
            estimiter = fminsearchout.iterations;
            
    end % END SWITCH
    
catch
    
    errcount=errcount+1;
    flag_optimok=0;
    
    try
        cd(outpath);
        diary errors;
        disp('Error in estima_loop');
        dispaj('Iteration Number ',ii);
        lasterr
        diary off;
        cd(cucd); 
    catch
        cd(backup_path);
        diary errors;
        disp('Error in estima_loop');
        dispaj('Iteration Number ',ii);
        lasterr
        diary off;
        cd(cucd);
    end % END TRY CATCH IN
    
end %END TRY CATCH OUT
cd(cucd); 

%cd(outpath); 
%save workspace; 
%cd(cucd); 

    
    
    time=toc/60;
    
    if flag_optimok==0
        continue
    end
    
    % =========================================================================
    % 8. OUTPUT
    % Output that must be produced by each code
    % -----------------------------------------
    % XESTMODE:  Estimated value
    % ESTIMITER: Number of Iterations
    % EXITFLAG : How routine ended
    % LPOSTMIN : Exit value of posterior (min)
    % =========================================================================
    parvecmode_temp = parvecmode; 
    parvecmode_temp(parposest) = xestmode;
    parestmode = xestmode;
    % =================================================================
    % 8.1. Evaluate Posterior and Likelihood at the MODE
    %      Write them to the temporary TAB matrix
    % =================================================================
    [lpostd, likel] = lmj_postmax(xestmode,parvecmode,parposest,funcmod,Y,trainvec, prpar, prss, solveopt, addsol, ssposest);
    
    logpostd_mat(ii, 2) = lpostd;
    loglikel_mat(ii, 2) = likel;
    logprior_mat(ii, 2) = lpostd - likel;
    
    % PARMODEMAT     Matrix of output vectors
    parmodemat(:, ii) = parvecmode_temp;
    
    
    optimat(:, ii) = [estimiter; exitflag; time]; 
    
%     optimat(1, ii) = estimiter;
%     optimat(2, ii) = exitflag;
%     optimat(3, ii) = time;
    
    % =================================================================
    % 8.2.
    % Print to screen
    % ===================================================================
    try
        disp(['Starting parameter ', num2str(ii), ': Estimated values are '])
        printcell([parnames(parposest), num2cprec(parestmode)])
        printcell([{'logPost'; 'logLht'}, num2cprec([lpostd; likel])])
    catch
        disp('errors in DISPLYAING RESULT')
    end
    
    % =================================================================
    % 8.3.
    % Write TAB_LOOP
    % ===================================================================
    
    tab_ii = [num2cprec(parmodemat(:,ii),15); num2str(logpostd_mat(ii,2)); ...
        num2str(loglikel_mat(ii,2));num2str(ii);num2str(estimiter);num2str(exitflag); ...
        num2str(time)]; 
    tab_loop{:, ii + 1} = tab_ii; 
    
%     tab_loop(1:numpar,ii+1)=num2cprec(parmodemat(:,ii),15);
%     tab_loop{numpar+2,ii+1}=num2str(logpostd_mat(ii,2));
%     tab_loop{numpar+3,ii+1}=num2str(loglikel_mat(ii,2));
%     tab_loop{numpar+4,ii+1}=num2str(ii);
%     tab_loop{numpar+5,ii+1}= num2str(estimiter); % Iteration count
%     tab_loop{numpar+6, ii+1}= num2str(exitflag); % Exit flag
%     tab_loop{numpar+7, ii+1}= num2str(time);
%     
    % =================================================================
    % 8.4.
    % Attemp to save output
    % ===================================================================
%     try
%         cd(outpath);
%         save mat_loop logprior_mat logpostd_mat logprior_mat ...
%             parmodemat optimat loglikel_mat
%         if ~isunix
%             xlswrite('tab_loop', tab_loop);
%         end
%         cd(cucd);
%     catch
%         cd(backup_path)
%         if ~isunix
%             xlswrite('tab_loop', tab_loop);
%         end
%         save mat_loop logprior_mat logpostd_mat logprior_mat ...
%             parmodemat optimat loglikel_mat
%         cd(cucd)
%         disp('Could not save the data')
%     end
    
    
    % =====================================================================
end

disp(' ');
disp('____________________________________');
disp('End all draws');

dispaj('Total time in minutes=',sum(optimat(3,:)));

if flag_diary==1
    diary off;
    copyfile(diaryfname,outpath);
    cd(cucd);
    eval(['delete ',diaryfname,'.txt']);
end

cd(outpath);
save mat_loop logprior_mat logpostd_mat logprior_mat ...
            parmodemat optimat loglikel_mat
if ~isunix
   xlswrite('tab_loop', tab_loop);
end
cd(cucd);


cd(outpath);
save workspace;
cd(cucd);

